Modelo de Solow

O modelo de Solow assume taxa de poupança consante, crescimento da população constante, e progresso tecnológico como exógeno. Um bem homogêneo é consumido e investido em capital $K$ tem técnologia de produção Cobb-Douglas \begin{equation} Y(t) = K(t)^\alpha ( A(t) L(t) )^{1-\alpha}, \quad 0 < \alpha < 1 \tag{1} \end{equation}

Suponha que o estoque de tecnologia $A(t)$ e o estoque de trabalho $L(t)$ ambos cresçam a taxas constantes: \begin{align} &\frac{dA(t)}{dt} = g, \quad &A(t) = A(t)\exp(gt); \ &\frac{dL(t)}{dt} = n, \quad &L(t) = L(t)\exp(nt). \end{align}

Podemos escrever (1) em termos de unidades efetivas de trabalho:

\begin{align} \frac{Y(t)}{A(t)L(t)} &= \frac{1}{(A(t)L(t))^\alpha (A(t)L(t))^{1-\alpha} } K(t)^\alpha (A(t)L(t))^{1-\alpha} \\ y(t) &= k(t)^\alpha \tag{2} \end{align}

Letras minúsculas denotam a variável em unidades de trabalho efetivo.

Definiremos a equação de movimento do capital efetivo por trabalhador (que descreve a trajetória da variação do estoque de capital, e, no nosso caso também do produto) como \begin{equation} \frac{dk(t)}{dt} = s y(t) - (n + \delta + g)k(t) \tag{3} \end{equation} \begin{equation} \frac{dk(t)}{dt} = s k(t)^\alpha - (n + \delta + g)k(t) \tag{4} \end{equation}

No estado-estacionário (steady-state) do modelo de Solow temos que $\frac{dk(t)}{dt}=0$. Note que isso implica $\frac{dy(t)}{dt}=0$ no Steady-State.

In [343]:
α = 0.33 # proporção do capital na renda
n = 0.02 # taxa de crescimento da população
δ = 0.01 # taxa de depreciação do capital físico
g = 0.04 # taxa de crescimento do produto
plot(plt_a, plt_b, grid=true, layout=@layout([a;b]))
Out[343]:

Igualando o lado direito da equação (4) a zero, temos que o valor para o qual o estoque de capital por unidade efetiva de trabalho converge no estado estacionário:

\begin{equation} k^* = \left( \frac{s}{n+g+\delta} \right)^{1/(1-\alpha)} \tag{5} \end{equation}

Tudo mais constante, $k^*$ é maior conforme a taxa de poupança cresce (mas cada vez menos!).

Agora vamos manipular equação (1) de modo a obter uma forma funcional que possa ser estimada por MQO Substiremos $k^*$ na função de produção (1), e dividiremos pelo estoque de trabalho $L(t)$, mas não pelo estoque de tecnologia $A(t)$ (por quê?) \begin{equation} Y(t) = K(t)^\alpha \left( A(t)L(t) \right)^{1-\alpha} \end{equation} \begin{equation} \frac{Y(t)}{L(t)} = \left( \frac{K(t)}{L(t)} \right)^\alpha \left( \frac{A(t)L(t)}{L(t)} \right)^{1-\alpha} \end{equation} \begin{equation} \frac{Y(t)}{L(t)} = \left( A(t)k(t) \right)^\alpha A(t)^{1-\alpha} \end{equation} \begin{equation} \frac{Y(t)}{L(t)} = (k^*)^\alpha A(t)^1 \end{equation} \begin{equation} \frac{Y(t)}{L(t)} = \left( \frac{s}{n+g+\delta} \right)^{\frac{\alpha}{1-\alpha}}A(0)\exp(gt) \end{equation} Tomando logaritmos: \begin{equation} \log \left( \frac{Y(t)}{L(t)} \right) = \log A(0) + gt + \frac{\alpha}{1-\alpha}\ln (s) - \frac{\alpha}{1-\alpha} \ln (n + g + \delta). \end{equation} Para contrastar a equação (6) com os dados temos que tratar $A(0)$. Não temos uma medida concreta da tecnologia inicial em cada país, então seguindo Mankiw et al. [1992] iremos fazer a seguinte suposição: \begin{equation} \log A(0) = \bar{a} + \epsilon_i, \end{equation} em $\bar{a}$ é uma constante e que $\epsilon_i$ é um choque país-específico. Finalmente, temos a equação (7) de Mankiw et al. [1992]: \begin{equation} \log \left( \frac{Y(t)}{L(t)} \right) = \bar{a} + gt + \frac{\alpha}{1-\alpha}\ln (s_i) - \frac{\alpha}{1-\alpha} \ln (n_i + g + \delta) + \epsilon_i. \tag{FE} \end{equation} Suponha $t=0$. Isto é: nosso "um" período $t=0$ vai de 1970 a 2014.

Algumas previsões do Modelo de Solow

(PA). A renda per capita de estado depende positivamente da taxa de poupança e negativamente da taxa de cresimento populacional. (previsões qualitativas)

(PB). Porque o modelo de Solow assume que os fatores de produção são remunerados de acordo com seus produtos marginais, e só temos dois fatores de produção; então sabemos que a magnitude da elasticidade-poupança da renda per-capita deve ser $\frac{\alpha}{1-\alpha}$, e que a magnitude da elasticidade-"n" da renda per-capita deve ser $-\frac{\alpha}{1-\alpha}$. (lembrando que $\alpha$ é a contribuição do capital para a renda) (previsão quantitativa).

Nosso objetivo é usar técnicas de regressão linear para testar (PA) e (PB).

In [344]:
df_vars = readtable("solow_vars.csv")[:,2:end];
head(df_vars)
Out[344]:
ISO3GDPPOPIOYGPOPSCHOOLdSCHOOLg
1AGO7967.857264324660.3390726996136360.030821547726231.200716850340910.008150951837462490.0419973074459723
2ALB10664.41940203680.1681053180.006320004224954652.413692980977270.01426360353410210.0329134983094553
3ARE64397.51251879370.3645012114772730.08132913018152512.150451969090910.01556446363477530.0502477150945378
4ARG20221.8293213690.161189991750.01319834042669372.534024027204550.007908883255217680.0503095429744743
5AUS43070.61938325960.2807338000454550.01362155172546293.355989862590910.003527761575123860.0322031572730742
6AUT47744.46084709540.2839526970681820.002822663578238162.995146577954550.005428562854949920.0306597636384636
In [345]:
plot(hst1,hst2,hst3,hst4)
Out[345]:
In [346]:
plot(lhst1,lhst2,lhst3,lhst4)
Out[346]:
In [347]:
sctr1
Out[347]:
In [348]:
sctr2
Out[348]:
In [349]:
sctr3
Out[349]:
In [351]:
df_vars_all = DataFrame()
df_vars_all[:y] = Array(log.(df_vars[:GDPPOP]))
df_vars_all[:ioy] = Array(log.(df_vars[:IOY]))
df_vars_all[:n] = Array(log.(df_vars[:GPOP]+0.05))
df_vars_all[:ioy_minus_n] = df_vars_all[:ioy]-df_vars_all[:n]
df_vars_all[:school] = log.(df_vars[:SCHOOL])
df_vars_all[:school_minus_n] = df_vars_all[:school]-df_vars_all[:n]
df_vars_all[:dschool] = log.(df_vars[:dSCHOOL]);

Forma Funcional 1: \begin{equation} \log \frac{Y_i}{L_i} = \beta_0 + \beta_1 \log s_i + \beta_2 \log(n_i + g + \delta) + \epsilon_i \end{equation}

In [491]:
# Estimação
y = Array(df_vars_all[:y])
X = Array([ones(y) df_vars_all[:ioy] df_vars_all[:n]]) # Primeira Especificação!
N = length(y) # tamanho amostral
p = size(X,2) # número de parâmetros estimados
β = inv(X'X)X'y
Out[491]:
3-element Array{Float64,1}:
  7.75598
  1.98036
 -1.79012
In [492]:
ϵ = y - X*β # resíduos
σ2_ϵ = ϵ'ϵ/(N-p) # variância dos resíduos
Out[492]:
0.7164298514227673
In [493]:
Σ2_β = σ2_ϵ*inv(X'X) # matriz de covariância do estimador β
Out[493]:
3×3 Array{Float64,2}:
 2.06777   0.143904   0.678709 
 0.143904  0.051206   0.0233547
 0.678709  0.0233547  0.237586 

Uma maneira de testar (PA) é um teste t individual para os parâmetros $\beta_1$ e $\beta_2$

$H_0$: $\beta_i = 0$

$H_1$: $\beta_i \neq 0$

Sabemos que $\frac{\beta_i}{SE(\beta_i)} \sim TDist_{(N-p)}$. Portanto podemos calcular as estatísticas de teste:

In [494]:
t_obs = β./sqrt.(diag(Σ2_β))
Out[494]:
3-element Array{Float64,1}:
  5.39369
  8.75154
 -3.67258
In [495]:
pvals = ccdf(TDist(N-p),abs.(t_obs))*2 #*2 porque Normal() é simétrica, e nosso teste é bicaudal.
Out[495]:
3-element Array{Float64,1}:
 4.7875e-7  
 6.16518e-14
 0.000391292
In [496]:
TSS = (y - mean(y))'*(y - mean(y)) #soma de quadrados totais
RSS = ϵ'ϵ # soma de quadrado dos resíduos
R2 = 1-RSS/TSS # coeficiente de determinação
adj_R2 = 1-(1-R2)*(N-1)/(N-p) # R² ajustado
df_out = DataFrame(); df_out[:coefs] = β; df_out[:t_obs] = t_obs; df_out[:pvals] = pvals; df_out[:R2]= R2; 
df_out[:adj_R2] = adj_R2
print(df_out)
3×5 DataFrames.DataFrame
│ Row │ coefs    │ t_obs    │ pvals       │ R2       │ adj_R2   │
├─────┼──────────┼──────────┼─────────────┼──────────┼──────────┤
│ 1   │ 7.75598  │ 5.39369  │ 4.7875e-7   │ 0.525551 │ 0.515868 │
│ 2   │ 1.98036  │ 8.75154  │ 6.16518e-14 │ 0.525551 │ 0.515868 │
│ 3   │ -1.79012 │ -3.67258 │ 0.000391292 │ 0.525551 │ 0.515868 │
In [497]:
reg1 = lm(@formula(y~ioy+n), df_vars_all)
Out[497]:
DataFrames.DataFrameRegressionModel{GLM.LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,Base.LinAlg.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}

Formula: y ~ 1 + ioy + n

Coefficients:
             Estimate Std.Error  t value Pr(>|t|)
(Intercept)   7.75598   1.43797  5.39369    <1e-6
ioy           1.98036  0.226287  8.75154   <1e-13
n            -1.79012  0.487428 -3.67258   0.0004

Forma Funcional 2: \begin{equation} \log \frac{Y_i}{L_i} = \gamma_0 + \gamma_1 [\log s_i - \log(n_i + g + \delta)] + \eta_i \end{equation}

In [498]:
reg2 = lm(@formula(y~ioy_minus_n), df_vars_all)
Out[498]:
DataFrames.DataFrameRegressionModel{GLM.LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,Base.LinAlg.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}

Formula: y ~ 1 + ioy_minus_n

Coefficients:
             Estimate Std.Error t value Pr(>|t|)
(Intercept)   7.28953  0.224388 32.4863   <1e-53
ioy_minus_n   1.93808  0.185267  10.461   <1e-16
In [500]:
# o "implied α" de Mankiw et al. [1992] na Tabela 1 vem da Forma Funcional 2 e da Forma Estrutural (FE).
β₁ = coef(reg2)[2]
implied_α = β₁/(1+β₁)
Out[500]:
0.6596421379901414

Vamos testar agora a previsão (PB): que os coeficientes de $log (s)$ e $log(n +g+\delta)$ têm mesma magnitude, mas sinais contrários. Para tal compararemos a acurácia dos dois modelos que já estimamos. Primeiro o modelo irrestrito (UR), no qual o parâmetros que multiplicam o log da taxa de poupança $log(s)$ e $log(n+g+\delta)$ podem ser diferentes. (nossa Forma Funcional 1!) \begin{equation} \log \frac{Y_i}{L_i} = \beta_0 + \beta_1 \log s_i + \beta_2 \log(n_i + g + \delta) + \epsilon_i \end{equation} Então o modelo restrito (R), no qual fixamos que os parâmetros têm mesma magnitude mas sinais contrários. (nossa Forma Funcional 2!) \begin{equation} \log \frac{Y_i}{L_i} = \gamma_0 + \gamma_1 [\log s_i - \log(n_i + g + \delta)] + \eta_i \end{equation}

Um Teste F:

$H_0$: $\beta_1$ = $\beta_2$ = $\gamma_1$.

$H_1$: caso contrário.

\begin{equation} F_{obs} = \frac{SSR_R - SSR_{UR}}{SSR_{UR}} \frac{N - q_{UR}}{q_R - q_{UR}} \sim F(q_{UR}-q_{R}, N-q_{UR}) \end{equation}

No nosso caso: $q_{UR} = 3$, $q_{R} = 2$, $N = 101$.

In [501]:
reg1 = lm(@formula(y ~ ioy + n), df_vars_all) # estimamos o modelo conforme a Forma Funcional 1
reg2 = lm(@formula(y ~ ioy_minus_n), df_vars_all) # estimamos o modelo conforme a Forma Funcional 2
resid1 = y - predict(reg1) # calculamos os resíduos da FF1
resid2 = y - predict(reg2) # calcularmos os resíduos da FF2
SSR_R = resid2'resid2 # Soma do quadrado dos resíduos da FF1
SSR_UR = resid1'resid1 # Soma do quadrado dos resíduos da FF2
q_UR = length(coef(reg1)) # Número de parâmetros na FF1
q_R = length(coef(reg2)) # Número de parâmetros na FF2
N = length(y) # Tamanho amostral 
F_obs = ((SSR_R-SSR_UR)/(q_R-q_UR))/(SSR_UR/(N-q_UR)) # Estatística do teste F que elaboramos
Out[501]:
-0.10787674618880172
In [502]:
pval_F = ccdf(FDist(q_UR-q_R, N-q_UR), F_obs) 
# pval = 1 indicando que a restrição de que β₁ = β₂ não piora significativamente o "fit" do modelo.
Out[502]:
1.0